1.1 Conditional Probability

To explain the concept of conditional probability we first look at some basic exercises in probability.

In the following questions, we are rolling two six-sided dice (denoted 2D6) with \(D_1\) and \(D_2\) denoting the result of the first and second die respectively, and we denote \(T\) as being the sum of the two - i.e.

\[ T = D_1 + D_2 \]

Question: We start simple - with no further detail, what is the probability of the total being 11?

Answer: We do not know anything about \(D_1\) or \(D_2\), so there are two outcomes that result in \(T = 11\):

\[(D_1 = 5, D_2 = 6) \text{ or } (D_1 = 6, D_2 = 5).\]

Thus, the probability of rolling 11 is

\[ P(T = 11) = \frac{2}{36} = 0.05556 \]

Question: What is the probability of getting 11 if the first dice is 5?

Answer: Out total is 11 if the second dice is 6, so there is only 1 possible outcome of those six possible rolls.

\[ P(T = 11 \mid D_1 = 5) = \frac{1}{6} = 0.1667 \]

1.2 Conditional Dependence

Now suppose we know the totals of each dice instead? Suppose we know that \(T = 9\) and we have the following distribution for \(D_1\):

\[ P(D_1) = \{0, \, 0, \, 0, \, 0, \, 0.5, \, 0.5 \} \] induces the following probability distribution for \(D_2\):

\[ P(D_2) = \{0, \, 0, \, 0.5, \, 0.5, \, 0, \, 0 \} \]

because \(D_1 = 5 \implies D_2 = 4\) and \(D_1 = 6 \implies D_2 = 3\)

Thus, while we conceptually think of \(D_1\) and \(D_2\) as being independent of one another, they become dependent if we have knowledge of \(T\).

In this case we say that \(D_1\) and \(D_2\) are conditionally dependent on one another given \(T\).

Now we define new variables, \(X_1\) and \(X_2\):

\[ X_1 = \begin{cases} 1 \text{ iff } T \text{ is even}\\ 0 \text{ otherwise} \end{cases} \; \; \; X_2 = \begin{cases} 1 \text{ iff } T >= 9\\ 0 \text{ otherwise} \end{cases} \]

Given no further information, we say that \(X_1\) and \(X_2\) are dependent, as knowing information about one variables affects our knowledge of the other.

However, what happens if we know something about \(T\)?

In this case, we say that \(X_1\) and \(X_2\) are independent given \(T\).

2.1 The Sprinkler Network

We start with a basic example of a Bayes Network: the sprinkler network.

yn <- c("Yes", "No")

cptable_lst <- list(
  cptable(~Rain,                        levels = yn, values = c(2, 8)),
  cptable(~Sprinkler + Rain,            levels = yn, values = c(1, 99, 4, 6)),
  cptable(~wetGrass + Rain + Sprinkler, levels = yn, values = c(99, 1, 8, 2, 9, 1, 0, 1))
)

sprinkler_cptlist <- cptable_lst %>% compileCPT()
sprinkler_grain   <- sprinkler_cptlist %>% grain()

sprinkler_grain %>% plot()

Two events can cause grass to be wet: Either the sprinkler is on or it is raining. Rain has a direct effect on the use of the sprinkler (namely that when it rains, the sprinkler is usually not turned on).

This can be modeled with a Bayesian network. The variables (R)ain, (S)prinkler, Wet(G)rass have two possible values: (y)es and (n)o.

We can factorise the joint probability mass function as

\[ p_{GSR}(g, s, r) = p_{G \mid SR}(g \mid s, r) p_{S \mid R}(s \mid r) p_R(r) \]

or overloading the notation a little:

\[ P(G, S, R) = P(G \mid S, R) \; P(S, R) = P(G \mid S, R) \; P(S \mid R) \; P(R) \]

This means we can construct the joint probability table by starting with the \(conditional probability tables\) (CPTs).

Create the 3 CPTs using the parray() function and the following conditional probabilities:

\[\begin{align*} P(R) &= 0.2 \\ P(S|R) &= 0.01 & P(S|\neg R) &= 0.4\\ P(G|S,R) &= 0.99 & P(G|S, \neg R) &= 0.9 & P(G|\neg S, R) &= 0.8 & P(G|\neg S, \neg R) &= 0 \end{align*}\]

First we want to construct this network using our various conditional probabilities:

yn <- c("yes", "no")

## P(R)
p_R <- parray(
  varNames = "Rain",
  levels = list(yn),
  values = c(0.2, 0.8)
)

## P(S|R)
p_S_R <- parray(
  varNames = c("Sprinkler", "Rain"),
  levels   = list(yn, yn), 
  values   = c(0.01, 0.99, 0.4, 0.6)
)

## P(G|S,R)
p_G_SR <- parray(
  varNames = c("GrassWet", "Sprinkler", "Rain"),
  levels   = list(yn, yn, yn),
  values   = c(0.99, 0.01, 0.8, 0.2, 0.9, 0.1, 0, 1)
)

ftable(p_G_SR, row.vars = "GrassWet")

##          Sprinkler  yes        no     
##          Rain       yes   no  yes   no
## GrassWet                              
## yes                0.99 0.90 0.80 0.00
## no                 0.01 0.10 0.20 1.00

We now combine these probabilities into a full joint distribution p_GSR:

p_GSR <- tabListMult(
  list(p_G_SR, p_S_R, p_R)
)

ftable(p_GSR, row.vars = "GrassWet")

##          Sprinkler     yes              no        
##          Rain          yes      no     yes      no
## GrassWet                                          
## yes                0.00198 0.28800 0.15840 0.00000
## no                 0.00002 0.03200 0.03960 0.48000

No suppose we know that the grass is wet. What is the probability it is raining?

p_RG  <- tabMarg(p_GSR, c("Rain", "GrassWet"))  ## P(R,G)
p_G   <- tabMarg(p_RG, "GrassWet")              ## P(G)
p_R_G <- tabDiv(p_RG, p_G)                      ## P(R|G)

Reading across this CPT, we see that if the grass is wet, there is about a 36% chance it is due to rain.

While the above methods of manipulating CPTs works, the package gRain provides us better functionality for answering questions such as this given observed evidence.

3.1 Querying the Network

By using the gRain package it allows us to use the in-built functionality to query the network.

To start with, we want to see the marginal distribution of alleles for the father.

genetic_family_grain %>%
  querygrain(nodes = "father")

## $father
## father
##   AA   AB   BB 
## 0.49 0.42 0.09

We can also look at the joint distribution of the mother and child:

genetic_family_grain %>%
  querygrain(
    nodes = c("child", "mother"),
    type  = "joint"
  )

##      mother
## child    AA    AB    BB
##    AA 0.343 0.147 0.000
##    AB 0.147 0.210 0.063
##    BB 0.000 0.063 0.027

We can pass this output to ftable() if we want:

genetic_family_grain %>%
  querygrain(
    nodes = c("child", "mother"),
    type  = "joint"
    ) %>%
  ftable(col.vars = "child")

##        child    AA    AB    BB
## mother                        
## AA           0.343 0.147 0.000
## AB           0.147 0.210 0.063
## BB           0.000 0.063 0.027

We can also produce the conditional distributions for each of the possible values for the father given values of the mother and child.

genetic_family_grain %>%
  querygrain(
    nodes = c("father", "child", "mother"),
    type  = "conditional"
    ) %>%
  ftable(col.vars = "father")

##              father   AA   AB   BB
## child mother                      
## AA    AA            0.70 0.30 0.00
##       AB            0.70 0.30 0.00
##       BB             NaN  NaN  NaN
## AB    AA            0.00 0.70 0.30
##       AB            0.49 0.42 0.09
##       BB            0.70 0.30 0.00
## BB    AA             NaN  NaN  NaN
##       AB            0.00 0.70 0.30
##       BB            0.00 0.70 0.30

3.2 Paternity Test

Now suppose we know that a child has genotype \(AB\), and the mother has genotype \(BB\). Given that a man has genotype \(AB\), what can we say about the chances of the man being the father of the child?

This is a more complicated question than it sounds, as we need to think about a few different ideas. We make use of the pEvidence() and setEvidence().

p_fmc <- genetic_family_grain %>%
  setEvidence(
    evidence = list(mother = "BB", child = "AB", father = "AB")
    )

p_fmc %>% print()

## Independence network: Compiled: TRUE Propagated: TRUE 
##   Nodes: chr [1:3] "child" "mother" "father"
##   Evidence:
##    nodes is.hard.evidence hard.state
## 1 mother             TRUE         BB
## 2  child             TRUE         AB
## 3 father             TRUE         AB
##   pEvidence: 0.018900

p_fmc %>% pEvidence()

## [1] 0.0189

So we now have a probability for finding this particular combination of alleles in a mother / father / child grouping.

So now we calculate the probability of any given man being \(AB\)?

p_f <- genetic_family_grain %>%
  setEvidence(evidence = list(father = "AB"))

p_f %>% print()

## Independence network: Compiled: TRUE Propagated: TRUE 
##   Nodes: chr [1:3] "child" "mother" "father"
##   Evidence:
##    nodes is.hard.evidence hard.state
## 1 father             TRUE         AB
##   pEvidence: 0.420000

p_f %>% pEvidence()

## [1] 0.42

Finally, we look at the probability of the child being \(AB\) with a mother \(BB\).

p_mc <- genetic_family_grain %>%
  setEvidence(
    evidence = list(mother = "BB", child = "AB")
    )

p_mc %>% print()

## Independence network: Compiled: TRUE Propagated: TRUE 
##   Nodes: chr [1:3] "child" "mother" "father"
##   Evidence:
##    nodes is.hard.evidence hard.state
## 1 mother             TRUE         BB
## 2  child             TRUE         AB
##   pEvidence: 0.063000

p_mc %>% pEvidence()

## [1] 0.063

So now we need to look at the various ratios of these probabilities to assess probabilities.

p_man <- pEvidence(p_fmc) / pEvidence(p_f)
print(p_man)

## [1] 0.045

p_father <- p_man / pEvidence(p_mc)
print(p_father)

## [1] 0.7142857

3.3 Extended Genetic Family

We can now extend this model to add extended families, such as uncle/aunts and grandparents.

c_mf <- parray(
  varNames = c("child", "mother", "father"),
  levels   = rep(list(genotypes), 3),
  values   = p_inheritance
  )

f_gmgf <- parray(
  varNames = c("father", "grandmother", "grandfather"),
  levels   = rep(list(genotypes), 3),
  values = p_inheritance
  )

u_gmgf <- parray(
  varNames = c("uncle", "grandmother", "grandfather"),
  levels   = rep(list(genotypes), 3),
  values   = p_inheritance
  )

m  <- parray("mother",      levels = list(genotypes), values = genotype_probs)
gf <- parray("grandfather", levels = list(genotypes), values = genotype_probs)
gm <- parray("grandmother", levels = list(genotypes), values = genotype_probs)


extended_family_grain <- list(c_mf, m, f_gmgf, u_gmgf, gm, gf) %>%
  compileCPT() %>%
  grain()

extended_family_grain %>% iplot()

Now suppose the man is dead and we do not have any genetic information for him, but we know his brother tests as \(AA\). How does this affect the evidence for potential paternity?

p_uncle <- extended_family_grain %>%
  setEvidence(evidence = list(mother = "BB", child = "AB", uncle = "AA"))

p_uncle %>% print()

## Independence network: Compiled: TRUE Propagated: TRUE 
##   Nodes: chr [1:6] "child" "mother" "father" "uncle" "grandmother" "grandfather"
##   Evidence:
##    nodes is.hard.evidence hard.state
## 1 mother             TRUE         BB
## 2  child             TRUE         AB
## 3  uncle             TRUE         AA
##   pEvidence: 0.037485

p_uncle %>% pEvidence()

## [1] 0.037485

We also want to show how we calculate similar values by using querygrain().

extended_family_grain %>%
  querygrain(
    nodes = c("child", "mother", "uncle"),
    type  = "joint"
    ) %>%
  ftable(col.vars = "uncle")

##              uncle       AA       AB       BB
## child mother                                 
## AA    AA           0.204085 0.123480 0.015435
##       AB           0.087465 0.052920 0.006615
##       BB           0.000000 0.000000 0.000000
## AB    AA           0.036015 0.082320 0.028665
##       AB           0.102900 0.088200 0.018900
##       BB           0.037485 0.022680 0.002835
## BB    AA           0.000000 0.000000 0.000000
##       AB           0.015435 0.035280 0.012285
##       BB           0.006615 0.015120 0.005265

4.1 Constructing the Bayesian Network

a    <- cptable(~ asia,
                values = c(1, 99),
                levels = yn)
t_a  <- cptable(~ tub | asia,
                values = c(5, 95, 1, 99),
                levels = yn)
s    <- cptable(~ smoke,
                values = c(5, 5),
                levels = yn)
l_s  <- cptable(~ lung | smoke,
                values = c(1, 9, 1, 99),
                levels = yn)
b_s  <- cptable(~ bronc | smoke,
                values = c(6, 4, 3, 7),
                levels = yn)
e_lt <- cptable(~ either | lung + tub,
                values = c(1, 0, 1, 0, 1, 0, 0, 1),
                levels = yn)
x_e  <- cptable(~ xray | either,
                values = c(98, 2, 5, 95),
                levels = yn)
d_be <- cptable(~ dysp | bronc + either,
                values = c(9, 1, 7, 3, 8, 2, 1, 9),
                levels = yn)

chestclinic_cpt_grain <- list(a, t_a, s, l_s, b_s, e_lt, x_e, d_be) %>%
  compileCPT() %>%
  grain()

chestclinic_cpt_grain %>% summary()

## Independence network: Compiled: TRUE Propagated: FALSE 
##  Nodes : chr [1:8] "asia" "tub" "smoke" "lung" "bronc" "either" "xray" "dysp"
##  Number of cliques:                 6 
##  Maximal clique size:               3 
##  Maximal state space in cliques:    8

We first should check some of the various conditional probabilities to ensure we have built the network correctly:

chestclinic_cpt_grain %>%
  querygrain(
    nodes = c("xray", "either"),
    type  = "conditional"
  )

##      either
## xray   yes   no
##   yes 0.98 0.05
##   no  0.02 0.95

chestclinic_cpt_grain %>%
  querygrain(
    nodes = c("dysp", "bronc", "either"),
    type  = "conditional"
    ) %>%
  ftable(row.vars = c("either", "bronc"))

##              dysp yes  no
## either bronc             
## yes    yes        0.9 0.1
##        no         0.7 0.3
## no     yes        0.8 0.2
##        no         0.1 0.9

4.2 Building Networks from Data

Now that we have our network, we can use data to construct the CPTs. Assuming no data is missing, we can use a tibble of data to set the probabilities.

data(chestSim500)
data(chestSim1000)
data(chestSim10000)


chestdata_500_tbl   <- chestSim500   %>% as_tibble()
chestdata_1000_tbl  <- chestSim1000  %>% as_tibble()
chestdata_10000_tbl <- chestSim10000 %>% as_tibble()

chestdata_500_tbl   %>% glimpse()

## Rows: 500
## Columns: 8
## $ asia   <fct> no, no, no, no, no, no, no, no, no, no, no, no, no, no, no, no…
## $ tub    <fct> no, no, no, no, no, no, no, no, no, no, no, no, no, no, no, no…
## $ smoke  <fct> no, no, yes, yes, yes, no, yes, no, no, no, no, no, yes, no, y…
## $ lung   <fct> no, no, no, no, no, no, no, no, no, no, no, no, yes, no, no, n…
## $ bronc  <fct> yes, yes, yes, yes, yes, no, yes, no, no, yes, no, no, yes, ye…
## $ either <fct> no, no, no, no, no, no, no, no, no, no, no, no, yes, no, no, n…
## $ xray   <fct> no, no, no, yes, no, no, no, no, no, no, no, no, yes, no, no, …
## $ dysp   <fct> no, yes, yes, no, yes, yes, yes, no, no, yes, no, no, yes, yes…

chestdata_1000_tbl  %>% glimpse()

## Rows: 1,000
## Columns: 8
## $ asia   <fct> no, no, no, no, no, no, no, no, no, no, no, no, no, no, no, no…
## $ tub    <fct> no, no, no, no, no, no, no, no, no, no, no, no, no, no, no, no…
## $ smoke  <fct> no, yes, yes, no, yes, yes, yes, no, yes, no, no, no, yes, yes…
## $ lung   <fct> no, no, no, no, no, yes, no, no, no, no, no, no, no, no, no, n…
## $ bronc  <fct> yes, yes, no, no, yes, yes, yes, no, no, no, no, no, yes, yes,…
## $ either <fct> no, no, no, no, no, yes, no, no, no, no, no, no, no, no, no, n…
## $ xray   <fct> no, no, no, no, no, yes, no, no, no, no, no, no, no, no, no, n…
## $ dysp   <fct> yes, yes, no, no, yes, yes, yes, no, no, no, no, no, yes, yes,…

chestdata_10000_tbl %>% glimpse()

## Rows: 10,000
## Columns: 8
## $ asia   <fct> no, no, no, no, no, no, no, no, no, yes, yes, no, no, no, no, …
## $ tub    <fct> no, no, no, no, no, no, no, no, no, no, no, no, no, no, no, no…
## $ smoke  <fct> no, yes, yes, yes, no, no, yes, yes, yes, no, no, yes, yes, ye…
## $ lung   <fct> no, no, no, no, no, no, no, no, no, no, no, no, no, no, no, no…
## $ bronc  <fct> no, yes, no, yes, no, yes, no, yes, yes, yes, yes, no, no, no,…
## $ either <fct> no, no, no, no, no, no, no, no, no, no, no, no, no, no, no, no…
## $ xray   <fct> no, no, no, no, no, no, no, no, no, no, no, no, no, no, no, no…
## $ dysp   <fct> no, yes, no, yes, no, yes, no, yes, no, yes, yes, no, no, no, …

We now can build multiple networks for each of the datasets we have.

chestsim_500_grain <- chestclinic_dag %>%
  grain(data = chestdata_500_tbl, smooth = 0.1) %>%
  compile(propagate = TRUE)

chestsim_1000_grain <- chestclinic_dag %>%
  grain(data = chestdata_1000_tbl, smooth = 0.1) %>%
  compile(propagate = TRUE)

chestsim_10000_grain <- chestclinic_dag %>%
  grain(data = chestdata_10000_tbl, smooth = 0.1) %>%
  compile(propagate = TRUE)

The unconditional probability of a patient having lung cancer should match the proportion of cases in the dataset.

chestsim_500_grain %>%
  querygrain(nodes = "lung", type = "joint")

## lung
##        yes         no 
## 0.04636229 0.95363771

chestdata_500_tbl %>%
  count(lung, name = "count") %>%
  mutate(prob = count / sum(count))

## # A tibble: 2 x 3
##   lung  count  prob
##   <fct> <int> <dbl>
## 1 yes      23 0.046
## 2 no      477 0.954

We do the same thing for the other networks we constructed.

First we look at the 1,000 row dataset.

chestsim_1000_grain %>%
  querygrain(nodes = "lung", type = "joint")

## lung
##        yes         no 
## 0.05317813 0.94682187

chestdata_1000_tbl %>%
  count(lung, name = "count") %>%
  mutate(prob = count / sum(count))

## # A tibble: 2 x 3
##   lung  count  prob
##   <fct> <int> <dbl>
## 1 yes      53 0.053
## 2 no      947 0.947

Now we look at the network based on 10,000 datapoints.

chestsim_10000_grain %>%
  querygrain(nodes = "lung", type = "joint")

## lung
##        yes         no 
## 0.05571778 0.94428222

chestdata_10000_tbl %>%
  count(lung, name = "count") %>%
  mutate(prob = count / sum(count))

## # A tibble: 2 x 3
##   lung  count   prob
##   <fct> <int>  <dbl>
## 1 yes     557 0.0557
## 2 no     9443 0.944

4.3 Using Evidence

We can use given evidence to adjust the conditional probabilities.

Suppose the individual has visited Asia and has dyspnoea, what is the conditional probability that the person has lung cancer?

chestsim_500_grain %>%
  querygrain(
    nodes    = "lung",
    evidence = list(asia = "yes", dysp = "yes"),
    type     = "marginal"
    )

## $lung
## lung
##        yes         no 
## 0.06414233 0.93585767

Given this evidence we see the probability of lung cancer has now increased to around 6.5%.

We repeat this exercise for the 1,000 network and see what probability it shows.

chestsim_1000_grain %>%
  querygrain(
    nodes    = "lung",
    evidence = list(asia = "yes", dysp = "yes"),
    type     = "marginal"
    )

## $lung
## lung
##       yes        no 
## 0.1030665 0.8969335

In this network, the probability is 10%.

chestsim_10000_grain %>%
  querygrain(
    nodes    = "lung",
    evidence = list(asia = "yes", dysp = "yes"),
    type     = "marginal"
    )

## $lung
## lung
##      yes       no 
## 0.103872 0.896128

As before, the probability is about 10%.

4.4 Marginal Proportions for Diseases

We also want to calculate the different marginal probabilities for the three conditions - Tuberculosis, Lung cancer and Dyspnoea - according to the different

chestclinic_cpt_grain %>%
  querygrain(
    nodes = c("lung", "bronc", "tub"),
    type  = "marginal"
    )

## $tub
## tub
##    yes     no 
## 0.0104 0.9896 
## 
## $lung
## lung
##   yes    no 
## 0.055 0.945 
## 
## $bronc
## bronc
##  yes   no 
## 0.45 0.55

chestsim_500_grain %>%
  querygrain(
    nodes = c("lung", "bronc", "tub"),
    type  = "marginal"
    )

## $tub
## tub
##        yes         no 
## 0.01432307 0.98567693 
## 
## $lung
## lung
##        yes         no 
## 0.04636229 0.95363771 
## 
## $bronc
## bronc
##       yes        no 
## 0.4540341 0.5459659

chestsim_1000_grain %>%
  querygrain(
    nodes = c("lung", "bronc", "tub"),
    type  = "marginal"
    )

## $tub
## tub
##         yes          no 
## 0.007197136 0.992802864 
## 
## $lung
## lung
##        yes         no 
## 0.05317813 0.94682187 
## 
## $bronc
## bronc
##       yes        no 
## 0.4360235 0.5639765

chestsim_10000_grain %>%
  querygrain(
    nodes = c("lung", "bronc", "tub"),
    type  = "marginal"
    )

## $tub
## tub
##        yes         no 
## 0.01091936 0.98908064 
## 
## $lung
## lung
##        yes         no 
## 0.05571778 0.94428222 
## 
## $bronc
## bronc
##       yes        no 
## 0.4422023 0.5577977

5.1 Scaling Up Bayesian Networks

All of the above approaches are example of the `Brute Force’ approach which is done by calculating the full joint distribution for the network \(p(V)\) as a multiple of the CPTs that comprise it.

\[ p(V) = p(a) \, p(t \mid a) \, p(s) \, p(l \mid s) \, p(b \mid s) \, p(e \mid t, l) \, p(d \mid e, b) \, p(x \mid e) \]

This gives \(p(V)\) represented by a table with \(2^8 = 256\) entries.

We then marginalise and condition as desired to calculate whatever probabilities we need.

The scaling of this approach is bad. A network with 80 variables, each with 10 values has a joint probability space of \(10^{80}\), approximately the count of atoms in the universe.

We are going to need a bigger boat…

So, we must find an approach that does not require the full joint distribution, instead focusing on the the low dimensional CPTs and send `messages’ between them.

To use a network it first needs to be compiled and then propagated.

Compilation of a network based on CPTs is first moralized — edges are added between the parents of each node, and then directed edges are replaced with undirected ones. It is then triangulated to form a triangulated graph.

The CPTs are transformed into clique potentials defined on the cliques of the chordal graph.

We see this process as below:

chestclinic_moralized    <- chestclinic_dag %>% moralize()
chestclinic_triangulated <- chestclinic_moralized %>% triangulate()

chestclinic_triangulated %>% plot()

Once we create the DAG, we view the triangulated data as a junction tree. The messages passed between connected nodes on the graph involve the common variables in the nodes, and propagating the information involves a double pass down and up the tree.

We now repeat the process from the above section, but now use the this approach to calculate our probabilities.

chestsim500_ev_grain <- chestsim_500_grain %>%
  setFinding(nodes = "asia", states = "yes", propagate = FALSE) %>%
  setFinding(nodes = "dysp", states = "yes", propagate = FALSE) %>%
  propagate()

We now use this network to perform the same calculations. As a reminder, this is calculating the marginal, joint and conditional probability for lung cancer and bronchitis given that the person recently visited Asia and exhibited shortness of breath (dyspnoea).

chestsim500_ev_grain %>%
  querygrain(
    nodes = c("lung", "bronc"),
    type  = "marginal"
    )

## $lung
## lung
##        yes         no 
## 0.06414233 0.93585767 
## 
## $bronc
## bronc
##       yes        no 
## 0.6761692 0.3238308

chestsim500_ev_grain %>%
  querygrain(
    nodes = c("lung", "bronc"),
    type  = "joint"
    )

##      bronc
## lung         yes         no
##   yes 0.03832565 0.02581668
##   no  0.63784353 0.29801414

chestsim500_ev_grain %>%
  querygrain(
    nodes = c("lung", "bronc"),
    type  = "conditional"
    )

##      bronc
## lung         yes         no
##   yes 0.05668056 0.07972275
##   no  0.94331944 0.92027725

5.2 Learning Network Structure

One final topic to think about is using the data to infer a model to explain it - so called model selection or structural learning.

Note that rather than using the data to determine the parameters of the model, we instead use the data to infer the structure of the model, that is, what the Bayesian network will look like.

For the chest clinic, we already know what the structure was, so we can compare the output of the various structural learning approaches to see how well it does.

Before we do that, we will plot the ‘true’ network:

chestclinic_dag %>% plot()

We now try a number of approaches to learn the structure - starting with the hill climbing algorithm implemented via hc().

chestSim500 %>%
  hc() %>%
  amat() %>%
  ggm::essentialGraph() %>%
  as("graphNEL") %>%
  plot()

We see that most of the network structure has been ‘learned’ from the data, apart from the relationship between “asia” and “tub” - presumably as this effect is small and thus weakly signalled within the data.

We repeat this exercise with chestSim10000 as this may allow us to pick that relationship out.

chestSim10000 %>%
  hc() %>%
  amat() %>%
  ggm::essentialGraph() %>%
  as("graphNEL") %>%
  plot()

We can try one larger set of data with 100,000 entries.

data(chestSim100000)

chestSim100000 %>%
  hc() %>%
  amat() %>%
  ggm::essentialGraph() %>%
  as("graphNEL") %>%
  plot()

6.1 Crop Analysis

We start with an example where we wish to analyse the crop yield of a particular plant. To do this, we have a number of things to consider:

1. The genetic and environmental potential of the plant
2. Vegatative mass
3. Harvest mass - the crop

For the first point, we construct two latent variables, \(G\) and \(E\) to measure the (G)enetic potential and (E)nvironmental potential of the plant.

We also want to measure various aspects of the (V)egatative mass, as these are indicators to the quality of the plant and hence affect the quality and quantity of the crop yield. Thus, we also have the variable \(V\).

The crop is a function of the vegatative mass and so is measured by the (N)umber of seeds and their mean (W)eight - giving us variables \(N\) and \(W\).

Finally, the (C)rop is determined by the number of seeds and the weight.

Converting all of the above into a Bayesian network, we can visualise it as follows:

(#fig:construct_crop_bn_dag)Plot of the Crop Network

If we were to specify a full joint probability distribution for \(p\) variables, this means we need \(p\) means, \(p\) variances and \(\frac{1}{2}p(p-1)\) elements of the correlation matrix and need to ensure that the correlation matrix is positive definite.

Thus, this problem scales as \(\mathcal{O}(N^2)\).

By using this graphical model, we only need to specify the various local, conditional distributions, simplifying the problem.

6.2 Constructing the Network

Before we add these distributions we want to construct the model. We can specify the above conditional factorisation of the model with the function model2network() function from the bnlearn package.

plant_crop_dagbnlearn <- model2network("[G][E][V|G:E][N|V][W|V][C|N:W]")

plant_crop_dagbnlearn %>% print()

## 
##   Random/Generated Bayesian network
## 
##   model:
##    [E][G][V|E:G][N|V][W|V][C|N:W] 
##   nodes:                                 6 
##   arcs:                                  6 
##     undirected arcs:                     0 
##     directed arcs:                       6 
##   average markov blanket size:           2.67 
##   average neighbourhood size:            2.00 
##   average branching factor:              1.00 
## 
##   generation algorithm:                  Empty

Suppose for example, we specify the following conditional distributions:

\[\begin{eqnarray*} G &\sim& \mathcal{N}(50, 10^2) \\ E &\sim& \mathcal{N}(50, 10^2) \\ V \mid G = g, E = e &\sim& \mathcal{N}(-10.35534 + 0.5 g + 0.70711 e, 5^2) \\ N \mid V = v &\sim& \mathcal{N}(45 + 0.1 v, 9.949874 ^ 2) \\ W \mid V = v &\sim& \mathcal{N}(15 + 0.7 v, 7.141428 ^ 2) \\ C \mid N = n, W = w &\sim& \mathcal{N}(0.3 n + 0.7 w, 6.25 ^ 2) \end{eqnarray*}\]

We now set up these conditional distributions in bnlearn. From inspection we see we are setting the parameters of a regression model at each node. As discussed, each node is distributed according to a Gaussian density, with conditional dependence manifesting as the conditional mean being a linear function of the dependent variables.

dist_E <- list(coef = c("(Intercept)" = 50), sd = 10)
dist_G <- list(coef = c("(Intercept)" = 50), sd = 10)
dist_V <- list(coef = c("(Intercept)" = -10.35534, E = 0.70711, G = 0.5), sd = 5)
dist_N <- list(coef = c("(Intercept)" = 45, V = 0.1), sd = 9.949874)
dist_W <- list(coef = c("(Intercept)" = 15, V = 0.7), sd = 7.141428)
dist_C <- list(coef = c("(Intercept)" = 0, N = 0.3, W = 0.7), sd = 6.25)

distrib_lst = list(E = dist_E, G = dist_G, V = dist_V, N = dist_N, W = dist_W, C = dist_C)


plant_crop_bnlearn <- custom.fit(plant_crop_dagbnlearn, dist = distrib_lst)

6.3 Inference on the Model

We now have constructed this ‘toy’ GBN and wish to focus on how we use this network and model to answer questions about the system - much like we did before in the discrete case.

Due to the use of Gaussian densities, we can use exact inference on this network. This has the benefit of calculating exact conditional probabilities, but imposes limits on the types of questions that can be ask.

As a result, we will focus on methods of approximate inference - that is, the use of Monte Carlo simulation to provide statistical estimates of the parameters of interest.

This gives us much more flexibility in terms of the types of questions we can ask, at the cost of producing answers that are only approximately correct.

That said, the error in this estimate is more often than not much smaller than any errors due to the choice of the model.

n_sim <- 1000000

plant_crop_sample_tbl <- rbn(n = n_sim, plant_crop_bnlearn) %>% as_tibble()

plant_crop_sample_tbl %>% glimpse()

## Rows: 1,000,000
## Columns: 6
## $ C <dbl> 44.80154, 57.24168, 65.23836, 62.43541, 53.70927, 64.48227, 50.8389…
## $ E <dbl> 57.81968, 41.31241, 48.04928, 55.16027, 53.30161, 58.30459, 64.3550…
## $ G <dbl> 49.37286, 63.04870, 72.86645, 36.11139, 47.21211, 48.66679, 56.3595…
## $ N <dbl> 48.00326, 71.30687, 60.71959, 48.65638, 46.91299, 45.58467, 57.6481…
## $ V <dbl> 52.10191, 49.71834, 66.85698, 55.27992, 56.06012, 56.21025, 58.1256…
## $ W <dbl> 55.09897, 57.57701, 62.78116, 67.73014, 56.41255, 60.74187, 42.8630…

ggplot(plant_crop_sample_tbl) +
  geom_histogram(aes(x = C), binwidth = 1) +
  scale_y_continuous(labels = label_comma()) +
  xlab("Value of C") +
  ylab("Simulation Count") +
  ggtitle("Marginal Distribution of C")

plant_crop_sample_tbl %>% filter(E == 50)

## # A tibble: 0 x 6
## # … with 6 variables: C <dbl>, E <dbl>, G <dbl>, N <dbl>, V <dbl>, W <dbl>

plant_crop_e50_5_tbl    <- plant_crop_sample_tbl %>% filter(abs(E - 50) < 5.0)
plant_crop_e50_1_tbl    <- plant_crop_sample_tbl %>% filter(abs(E - 50) < 1.0)
plant_crop_e50_01_tbl   <- plant_crop_sample_tbl %>% filter(abs(E - 50) < 0.1)
plant_crop_e50_001_tbl  <- plant_crop_sample_tbl %>% filter(abs(E - 50) < 0.01)

plant_crop_e50_tbl <- list(
    `Width 5.0`  = plant_crop_e50_5_tbl,
    `Width 1.0`  = plant_crop_e50_1_tbl,
    `Width 0.1`  = plant_crop_e50_01_tbl,
    `Width 0.01` = plant_crop_e50_001_tbl
    ) %>%
  bind_rows(.id = "width_label")

plot_tbl <- plant_crop_e50_tbl %>%
  group_by(width_label) %>%
  summarise(
    .groups = "drop",
    
    mean_value = mean(C),
    
    p10 = quantile(C, 0.10),
    p25 = quantile(C, 0.25),
    p50 = quantile(C, 0.50),
    p75 = quantile(C, 0.75),
    p90 = quantile(C, 0.90)
    )

ggplot(plot_tbl) +
  geom_errorbar(aes(x = width_label, ymin = p25, ymax = p75),
                width = 0, size = 3, colour = "red") +
  geom_errorbar(aes(x = width_label, ymin = p10, ymax = p90),
                width = 0, size = 1, colour = "black") +
  geom_point(aes(x = width_label, y = p50),        colour = "red") +
  geom_point(aes(x = width_label, y = mean_value), colour = "black") +
  expand_limits(y = 0) +
  xlab("Interval Width") +
  ylab("Value of C") +
  ggtitle("Conditional Values of C")

ggplot(plant_crop_e50_tbl) +
  geom_histogram(aes(x = C), binwidth = 1) +
  facet_wrap(vars(width_label), scales = "free_y", ncol = 2) +
  xlab("C") +
  ylab("Frequency Count") +
  ggtitle("Facetted Histogram Plot")

ggplot(plant_crop_e50_tbl) +
  geom_line(aes(x = C, colour = width_label), stat = "density") +
  xlab("C") +
  ylab("Density") +
  ggtitle("Comparison Plot of Estimated Density Functions")

plant_crop_e70_tbl <- list(
    `Width 5.0`  = plant_crop_sample_tbl %>% filter(abs(E - 70) < 5.0),
    `Width 1.0`  = plant_crop_sample_tbl %>% filter(abs(E - 70) < 1.0),
    `Width 0.1`  = plant_crop_sample_tbl %>% filter(abs(E - 70) < 0.1),
    `Width 0.01` = plant_crop_sample_tbl %>% filter(abs(E - 70) < 0.01)
    ) %>%
  bind_rows(.id = "width_label")

plant_crop_e70_tbl %>% glimpse()

## Rows: 72,093
## Columns: 7
## $ width_label <chr> "Width 5.0", "Width 5.0", "Width 5.0", "Width 5.0", "Widt…
## $ C           <dbl> 55.16532, 65.78010, 53.84612, 52.69393, 56.45567, 61.6004…
## $ E           <dbl> 68.56438, 67.05995, 71.12048, 66.88804, 71.49582, 69.3356…
## $ G           <dbl> 54.32818, 56.42899, 52.84883, 55.80996, 41.14224, 41.3920…
## $ N           <dbl> 52.35925, 57.40292, 50.67462, 60.10386, 59.95746, 43.6093…
## $ V           <dbl> 67.22984, 67.47149, 68.79282, 57.69178, 62.78158, 54.2236…
## $ W           <dbl> 61.06274, 56.41881, 66.60216, 59.49617, 52.83071, 63.2764…

ggplot(plant_crop_e70_tbl) +
  geom_line(aes(x = C, colour = width_label), stat = "density") +
  xlab("C") +
  ylab("Density") +
  ggtitle("Comparison Plot of Estimated Density Functions")

plant_crop_goodcrop_tbl <- plant_crop_sample_tbl %>%
  filter(C > 80)

ggplot(plant_crop_goodcrop_tbl) +
  geom_point(aes(x = E, y = G)) +
  xlab("E") +
  ylab("G") +
  ggtitle("Conditional Joint Distribution of E and G")

plant_crop_importance_tbl <- cpdist(
    plant_crop_bnlearn,
    nodes    = c("E", "G"),
    evidence = (C > 80),
    n = 1000000
    ) %>%
  as_tibble()

plant_crop_importance_tbl %>% glimpse()

## Rows: 1,342
## Columns: 2
## $ E <dbl> 47.46442, 68.75327, 67.12545, 60.86390, 48.37759, 59.48258, 67.9352…
## $ G <dbl> 60.10531, 60.65128, 48.64428, 53.58544, 55.32707, 67.23726, 78.7154…

ggplot(plant_crop_importance_tbl) +
  geom_point(aes(x = E, y = G)) +
  xlab("Environmental Potential") +
  ylab("Genetic Potential") +
  ggtitle("Plot of High Crop Yield Environmental Variables")

comparison_tbl <- list(
    simulation = plant_crop_goodcrop_tbl,
    importance = plant_crop_importance_tbl
    ) %>%
  bind_rows(.id = "source")

ggplot(comparison_tbl) +
  geom_point(aes(x = E, y = G, colour = source)) +
  xlab("Environmental Potential") +
  ylab("Genetic Potential") +
  ggtitle("Comparison Plot of (E,G) Where C > 80")